R Bootcamp day 3
Statistical Computing Track
trojan R logo

George G. Vega Yon


University of Southern California
Department of Preventive Medicine

August 21th, 2019

Before we start

  1. You need to have Rcpp installed in your system:

  2. You need to have a compiler

    • Windows: You can download Rtools from here.

    • MacOS: It is a bit complicated… Here are some options:

      • CRAN’s manual to get the clang, clang++, and gfortran compilers here.

      • A great guide by the coatless professor here

And that’s it!

R is great, but…

R is great, but… (cont’d)

Enter Rcpp

Dirk Romain

Dirk Eddelbuettel and Romain François, lead authors of this R package

From the package description:

The ‘Rcpp’ package provides R functions as well as C++ classes which offer a seamless integration of R and C++

Why bother?

Example 1: Looping over a vector

#include<Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector add1(NumericVector x) {
  NumericVector ans(x.size());
  for (int i = 0; i < x.size(); ++i)
    ans[i] = x[i] + 1;
  return ans;
}
add1(1:10)
##  [1]  2  3  4  5  6  7  8  9 10 11

Example 1: Looping over a vector (vers 2)

Make it more sweet by adding some “sugar” (the Rcpp kind)

#include<Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector add1Cpp(NumericVector x) {
  return x + 1;
}
add1Cpp(1:10)
##  [1]  2  3  4  5  6  7  8  9 10 11

How much fast?

Compared to this:

add1R <- function(x) {
  for (i in 1:length(x))
    x[i] <- x[i] + 1
  x
}

microbenchmark::microbenchmark(add1R(1:1000), add1Cpp(1:1000), unit = "relative")
## Unit: relative
##             expr      min       lq    mean   median       uq      max
##    add1R(1:1000) 22.62702 23.05899 10.1886 21.60072 17.68509 4.660335
##  add1Cpp(1:1000)  1.00000  1.00000  1.0000  1.00000  1.00000 1.000000
##  neval cld
##    100   b
##    100  a

Main differences between R and C++

  1. One is compiled and the other interpreted

  2. Indexing objects: In C++ the indices range from 0 to (n - 1), whereas in R is from 1 to n.

  3. All expressions end with a ; (this is optional in R).

  4. In C++ object need to be declared, in R not (dynamic).

C++/Rcpp fundamentals: Types

Besides of the C-like types (double, int, char, and bool), we can use the following types of objects with Rcpp:

C++/Rcpp fundamentals: Parts of “an Rcpp program”

Line by line, we see the following:

  1. The #include<Rcpp.h> is similar to library(...) in R, it brings in all that we need to write C++ code for Rcpp.

  2. using namespace Rcpp is somewhat similar to detach(...). This simplifies syntax. If we don’t include this, all calls to Rcpp members need to be explicit, e.g. instead of typing NumericVector, we would need to type Rcpp::NumericVector

  3. The // starts a comment in C++, in this case, // [[Rcpp::export]] comment is a flag that Rcpp uses to “export” this C++ function to R.

  4. Is the first part of the function definition. We are creating a function that returns a NumericVector, is called add1, has a single input element named x that is also a NumericVector.

C++/Rcpp fundamentals: Parts of “an Rcpp program” (cont’d)

  1. Here we are declaring an object called ans which is a NumericVector with initial size equal to the size of x. Notice that .size() is called a “member function” of the x object, which is of class NumericVector.

  2. We are declaring a for-loop. This has three parts:

    1. int i = 0 We declare the variable i, an integer, and initialize it at 0.

    2. i < x.size() This loop will end when i’s value is at or above the length of x.

    3. ++i At each iteration, i will increment in one unit.

  3. ans[i] = x[i] + 1 set the i-th element of ans equal to the i-th element of x plus 1.

  4. return ans exists the function returning the vector ans.

C++/Rcpp fundamentals (cont’d 2)

Now, where to execute/run this?

For now, let’s just use the first option.

Example running .cpp file

Imagine that we have the following file named norm.cpp

#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
double normRcpp(NumericVector x) {
  
  return sqrt(sum(pow(x, 2.0)));
  
}

We can compile and obtain this function using this line Rcpp::sourceCpp("norm.cpp"). Once compiled, a function called normRcpp will be available in the current R session.

Now, get ready for some Rcpp action!

Problem 1: Adding vectors

  1. Using what you have just learned about Rcpp, write a function to add two vectors of the same length. Use the following template
```cpp
#include <Rcpp.h>

using namespace Rcpp;

// [[Rcpp::export]]
NumericVector add_vectors([declare vector 1], [declare vector 2]) {
  
  ... magick ...
  
  return [something];
}
```
    
  1. Now, we have to check for lengths. Use the stop function to see how can you return with an error if the two vectors don’t match size. Add the following lines in your code
```cpp
if ([some condition])
  stop("an arbitrary error message :)");
```

Problem 2: Fibonacci series

Fibonacci Spiral

$$
F(n) = \left\{\begin{array}{ll}
n, & \mbox{ if }n \leq 1\\
F(n - 1) + F(n - 2), & \mbox{otherwise}
\end{array}\right.
$$
```r
fibR <- function(n) {
  if (n <= 1)
    return(n)
  fibR(n - 1) + fibR(n - 2)
}

# Is it working?
c(fibR(0), fibR(1), fibR(2), fibR(3), fibR(4), fibR(5), fibR(6))
```

```
## [1] 0 1 1 2 3 5 8
```

Now, let’s translate this code into Rcpp and see how much speed boost we get!

Problem 2: Fibonacci series (solution)

#include <Rcpp.h>

// [[Rcpp::export]]
int fibCpp(int n) {
  if (n <= 1)
    return n;
  
  return fibCpp(n - 1) + fibCpp(n - 2);
  
}
microbenchmark::microbenchmark(
  fibR(20),
  fibCpp(20), unit="relative"
)
## Unit: relative
##        expr      min       lq     mean   median       uq      max neval
##    fibR(20) 257.8901 256.8613 182.5437 230.5804 205.0108 12.12437   100
##  fibCpp(20)   1.0000   1.0000   1.0000   1.0000   1.0000  1.00000   100
##  cld
##    b
##   a

Thank you!

R Bootcamp day 3
Statistical Computing Track
trojan R logo

website: ggvy.cl

twitter: @gvegayon

github: gvegayon

Session info

devtools::session_info()
## ─ Session info ──────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       Ubuntu 16.04.6 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language en_US                       
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Los_Angeles         
##  date     2019-08-21                  
## 
## ─ Packages ──────────────────────────────────────────────────────────────
##  package        * version    date       lib source                        
##  assertthat       0.2.1      2019-03-21 [1] CRAN (R 3.5.3)                
##  backports        1.1.4      2019-04-10 [1] CRAN (R 3.6.0)                
##  callr            3.3.1      2019-07-18 [1] CRAN (R 3.6.1)                
##  cli              1.1.0      2019-03-19 [1] CRAN (R 3.5.3)                
##  codetools        0.2-16     2018-12-24 [3] CRAN (R 3.5.2)                
##  crayon           1.3.4      2017-09-16 [1] CRAN (R 3.4.1)                
##  desc             1.2.0      2018-05-01 [1] CRAN (R 3.4.4)                
##  devtools         2.1.0      2019-07-06 [1] CRAN (R 3.6.0)                
##  digest           0.6.20     2019-07-04 [1] CRAN (R 3.6.0)                
##  evaluate         0.14       2019-05-28 [1] CRAN (R 3.6.0)                
##  fs               1.3.1      2019-05-06 [1] CRAN (R 3.6.0)                
##  glue             1.3.1      2019-03-12 [1] CRAN (R 3.5.3)                
##  htmltools        0.3.6      2017-04-28 [1] CRAN (R 3.5.1)                
##  knitr            1.24       2019-08-08 [1] CRAN (R 3.6.1)                
##  lattice          0.20-38    2018-11-04 [1] CRAN (R 3.5.1)                
##  magrittr         1.5        2014-11-22 [1] CRAN (R 3.5.1)                
##  MASS             7.3-51.4   2019-04-26 [3] CRAN (R 3.6.0)                
##  Matrix           1.2-17     2019-03-22 [3] CRAN (R 3.5.3)                
##  memoise          1.1.0      2017-04-21 [1] CRAN (R 3.4.0)                
##  microbenchmark   1.4-6      2018-10-18 [1] CRAN (R 3.5.1)                
##  multcomp         1.4-10     2019-03-05 [1] CRAN (R 3.5.3)                
##  mvtnorm          1.0-11     2019-06-19 [1] CRAN (R 3.6.0)                
##  pkgbuild         1.0.4      2019-08-05 [1] CRAN (R 3.6.1)                
##  pkgload          1.0.2      2018-10-29 [1] CRAN (R 3.5.1)                
##  prettyunits      1.0.2      2015-07-13 [1] CRAN (R 3.4.1)                
##  processx         3.4.1      2019-07-18 [1] CRAN (R 3.6.1)                
##  ps               1.3.0      2018-12-21 [1] CRAN (R 3.5.3)                
##  R6               2.4.0      2019-02-14 [1] CRAN (R 3.5.3)                
##  Rcpp             1.0.2      2019-07-25 [1] CRAN (R 3.6.1)                
##  remotes          2.1.0      2019-06-24 [1] CRAN (R 3.6.0)                
##  rlang            0.4.0      2019-06-25 [1] CRAN (R 3.6.0)                
##  rmarkdown        1.14       2019-07-12 [1] CRAN (R 3.6.1)                
##  rprojroot        1.3-2      2018-01-03 [1] CRAN (R 3.4.3)                
##  sandwich         2.5-1      2019-04-06 [1] CRAN (R 3.6.0)                
##  sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 3.5.1)                
##  stringi          1.4.3      2019-03-12 [1] CRAN (R 3.5.3)                
##  stringr          1.4.0      2019-02-10 [1] CRAN (R 3.5.3)                
##  survival         2.44-1.1   2019-04-01 [3] CRAN (R 3.5.3)                
##  testthat         2.2.1      2019-07-25 [1] CRAN (R 3.6.1)                
##  TH.data          1.0-10     2019-01-21 [1] CRAN (R 3.5.3)                
##  usethis          1.5.1.9000 2019-07-12 [1] Github (r-lib/usethis@c3ad943)
##  withr            2.1.2      2018-03-15 [1] CRAN (R 3.4.3)                
##  xfun             0.9        2019-08-21 [1] CRAN (R 3.6.1)                
##  yaml             2.2.0      2018-07-25 [1] CRAN (R 3.5.1)                
##  zoo              1.8-6      2019-05-28 [1] CRAN (R 3.6.0)                
## 
## [1] /usr/local/lib/R/site-library
## [2] /usr/lib/R/site-library
## [3] /usr/lib/R/library